arXiv: 1508.03137v2 [physics.geo-ph] 16 Oct 2015 


Torsion of a cylinder of partially molten rock with a spherical 
inclusion: theory and simulation 

Laura Alisic* Sander Rhebergen’^ John F. Rudge^ Richard F. Katz^ 

Garth N. Wells^ 


Abstract 

The processes that are involved in migration and extraction of melt from the mantle 
are not yet fully understood. Gaining a better understanding of material properties of 
partially molten rock could help shed light on the behavior of melt on larger scales in the 
mantle. In this study, we simulate three-dimensional torsional deformation of a partially 
molten rock that contains a rigid, spherical inclusion. We compare the computed porosity 
patterns to those found in recent laboratory experiments. The laboratory experiments show 
emergence of melt-rich bands throughout the rock sample, and pressure shadows around the 
inclusion. The numerical model displays similar melt-rich bands only for a small bulk-to- 
shear-viscosity ratio (five or less). The results are consistent with earlier two-dimensional 
numerical simulations; however, we show that it is easier to form melt-rich bands in three 
dimensions compared to two. The addition of strain-rate dependence of the viscosity causes 
a distinct change in the shape of pressure shadows around the inclusion. This change in 
shape presents an opportunity for experimentalists to identify the strain-rate dependence 
and therefore the dominant deformation mechanism in torsion experiments with inclusions. 


1 Introduction 


The transport of melt in the mantle plays an important role in the dynamics and chemical evo¬ 
lution of both the mantle and the crust. Although the equations that describe the conservation 


of mass, momentum, and energy of partially molten rock are well established [McKenzie, 1984 


Bercovici and Ricard, 2003), the appropriate constitutive relations remain uncertain. This means 


that the dynamics of melt segregation and transport present significant unanswered questions. 

One means of addressing questions on the dynamics of melt segregation and transport is 
by comparison of simulations with laboratory experiments on partially molten rocks subjected 
to forced deformation. A recent experimental study with significant potential in this regard is 


reported by Qi et al (2013). Following on the torsional deformation experiments of King et al. 


(2010), Qi et al (2013) modified the basic experiment by including rigid, spherical beads within 
the partially molten rock that is undergoing deformation. They find that pressure shadows 
around the bead are expressed as variations in melt fraction there. Furthermore, they find 
that melt-rich bands, also observed in experiments without beads (e.g. Holtzman et a/., 2003), 
emerge and tend to connect with the large-porosity lobes of the pressure shadow. Previous 
analysis of pressure shadows [McKenzie and Holness, 2000| Rudge, 2014) and their interaction 
with banding instabilities [Alisic et a/., 2014) suggests that the observed relationship between 
these two modes of compaction could constrain the bulk viscosity of the two-phase system. 
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In working towards a better understanding of compaction in a two-phase system, we pose the 
following questions about the viscosity of the two-phase system that remain unresolved. What 


is the ratio of the bulk viscosity to the shear viscosity at small reference porosity (Simpson 


et a/., 2010)? How do the bulk and shear viscosities vary with porosity (e.g. Kelemen et al 
1997; Mei et a/., 2002; Takei and Holtzman, 2009)7 Is the rheology non-Newtonian and, if so. 


does this help to explain the patterns observed in experiments? And, more broadly, is a solely 
viscous rheology sufficient to capture the dynamics? These are long-term questions that we 
address. However, we find that on the basis of the comparison between experiments and theory 
considered here, we cannot answer these questions definitively, and we present a discussion of 
this shortcoming. 

In an earlier paper, we developed two-dimensional models of two-phase flow around a cylin¬ 
drical inclusion to study the same experimental system (Alisic et al, 2014). Here we build on 


those results by expanding the numerical simulations to three dimensions. This allows us to 
capture the three-dimensional scaling of compaction around a sphere, which differs from the 
two-dimensional scaling around a cylinder (Rudge, 2014). Moreover, the simulations presented 


here provide a more realistic comparison to the results of laboratory experiments ( Qi et al 
2013). These new simulations with ^ 7 x 10^ degrees of freedom would be impossible without 


an advanced, new preconditioning method for the equations of magma dynamics that has been 


recently developed (Rhebergen et a/., 2015). 


We begin this manuscript with a description of the equations governing deformation and 
compaction of partially molten rock, after which we summarize the domain geometry, boundary 
conditions, and discretization used in the numerical simulations. Analytical solutions for certain 
limiting cases are provided in Appendix we use these to benchmark the simulation code. 
The first set of results in Section |3.1| pertains to simulations with a uniform initial porosity 
that allow us to focus on compaction around a spherical inclusion in three dimensions, with 
Newtonian and non-Newtonian rheology. The simulations in Section 3.2 focus on problems with 
a random initial porosity field, where we investigate the interaction between pressure shadows 
around the inclusion and melt-rich bands developing throughout the domain. The results are 
followed by a discussion in Section]^ after which conclusions are drawn. 


2 The model 


2.1 Governing equations 


The compaction of partially molten rock and the transport of melt can be described by governing 
equations for two-phase flow, formulated here following MeKenzie (1984). In dimensionless form 
(see Appendix for the nondimensionalization): 


dd) , . 

-V ■ u. + V ■ 

-V Us - {RC,)~^Pc 
— V • T -h Vpf + Vpc 


= 0 , 

( 1 ) 

= 0 , 

( 2 ) 

= 0 , 

( 3 ) 

= 0 , 

( 4 ) 


where t denotes time, (j) is porosity, is the solid (matrix) velocity, pf and Pc are the magma 
and compaction pressure, respectively, and r is the deviatoric stress in the solid. Constitutive 
properties, discussed further below, appear as for the permeability and for the bulk vis¬ 
cosity. The bulk-to-shear-viscosity ratio in the reference state is defined as i? = Cref/'^ref? where 
is a reference bulk viscosity and r/ref is a reference shear viscosity for the two-phase mixture. 
Finally, D = 6/H where 6 is the compaction length and H is the height of the domain. The 
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compaction length is given by (McKenzie, 1984): 


5 = 
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(5) 


where Kref is the permeability in the reference state, and /r/ is the magma viscosity. In this 
study we assume a compaction length that is much larger than the domain size {D = 100). The 
deviatoric stress tensor r is 


T = r][ Vus + (Vus 




Us 


( 6 ) 


2e 


where r] is the shear viscosity and e is the deviatoric strain-rate tensor. 

The above model assumes that no melting or solidification takes place, buoyancy forces are 
negligible, and that the fluid and solid phases have densities that are constant (but different 
from each other); these assumptions are appropriate for the motivating laboratory experiments 


(e.g. Holtzman et al , 2003), though a model for the Earth’s mantle clearly must be more general. 


The unknown fields in the model are 0, u^, pj, and Pc? which must satisfy equations 0-(i> 
subject to the boundary conditions described below. 


2.2 Rheology and permeability 

Closure conditions for the governing equations are prescribed as 

K^, = (^^) , 7] = C = (^) 


( 7 ) 


where n and m depend on the melt geometry considered. We take n — 2 and m = 1, assuming 
a tubular melt geometry. In these definitions, 00 is the reference porosity, a is a, constant 
representing the sensitivity of matrix shear viscosity to porosity, ^ is the second invariant of the 
deviatoric strain-rate tensor, 

/I 

^ = ( : e ) , (8) 


and q is related to the power-law exponent n by 
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A power-law exponent n = 1 gives the limit of Newtonian rheology. 

In this study we focus on the effects of the porosity sensitivity cr, the reference bulk-to- 
shear-viscosity ratio i?, and the power-law exponent n on compaction patterns around and 
away from an inclusion. Laboratory experiments indicate that a is around 26 for a Newtonian 


et al. 


1997; Mei 


rheology (diffusion creep) and around 31/n for dislocation creep (Kelemen et al 

2002[). Sever al previous modelling studies have used a value cr = 28 (e.g. Alisie et al. 


\2014:) ^Katz et al (2006)). The bulk-to-shear-viscosity ratio, however, is poorly constrained. 


Theoretical and experimental studies place R between order one, independent of the reference 
porosity ( Takei and Holtzman\ 2009), and ^ 20 for a reference porosity 0o = 0.05 on the 


basis of the expected relation R oc 0“^ (Bereoviei and Rieard, 2003; Simpson et a/., 2010). In 
the simulations presented here, we vary a between zero and 50, and R between 5/3 and 20. In 
simulations where we study the effect of strain-rate dependence of the shear and bulk viscosities, 
the power-law exponent n has values between 1 and 6. The deformation mechanism of diffusion 
creep corresponds to n = 1, resulting in Newtonian viscosities. A larger exponent of around 3 
to 4 is relevant for dislocation creep. Katz et al. (2006) showed that increasing the power-law 


exponent is one possible mechanism for reproducing the shallow angle of melt-rich bands as 
observed in laboratory experiments. We therefore include simulations with n up to six. 
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Figure 1: The geometry of the model domain: a cylinder of height H — 1 and radius one (only 
a thin cut-out is shown in light gray), with a rigid spherical inclusion of radius 0 . 1 . On a two- 
dimensional slice through the cylinder and inclusion at x = the instantaneous compaction 
rate (V • u^) is plotted for a simulation with bulk-to-shear-viscosity ratio i? = 20, and power-law 
exponent n = 1 at time t = 0. The arrows at the top, bottom and side of the cylinder indicate 
the prescribed solid velocity on the cylinder boundaries. 

2.3 Domain of interest and boundary and initial conditions 

We compute solutions to equations 0-0 in a cylindrical domain C of height H — 1 
and radius 1 , where < 1 and 0 < z < 1 . A rigid spherical inclusion is centered at 

ro = (1/2,0,1/2) and has a radius of 0.1, as shown in Figure The inclusion is modeled as a 
spherical hole in the domain that cannot deform. The boundary conditions are defined on the 
boundary dfl as: 


KfpVpf • n = 0 on dfl, 
U 5 = w on dfl. 


( 10 ) 

( 11 ) 


where the boundaries are taken to be impermeable (equation ®), and w is a prescribed solid 
velocity. The cylinder is placed under torsion. This torsion is enforced by Dirichlet boundary 
conditions of the form 0 on the top and bottom (z = 0 and 2: = 1 ), and side boundaries of 
the cylinder (x^ -h = 1), such that on the outside of the cylinder w = Ucyp 


Ucyl — 


(-y (^ - 5) ’ ® (^ - 5) ’ 0) • 


( 12 ) 


Formally, the boundary conditions on the rigid inclusion are conditions of no net force and no net 
torque (see Alisic et aL ( 2014[ Appendix B)). Here, to simplify the construction of the numerical 
model, we instead apply a Dirichlet boundary condition of the form 0 that approximates the 
zero net force and torque conditions. In a uniform medium, a rigid sphere placed at the mid¬ 
plane of a torsion field should not translate, but should rotate with an angular velocity equal 
to half the vorticity of the imposed torsion field (see Section B.l). Thus we use a Dirichlet 
condition on the boundary of the inclusion with w = Ugphere? 


^sphere X (r Tq) , 


(13) 


where r = {x, y, z) is the position vector, and ro is the center of the sphere. The angular velocity 
n is given by 

^ = -V X Ucyl I 1 , 0 , 0 ^ . ( 14 ) 
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In our Cartesian coordinate system and for the position of the inclusion, (13) can thus be written 


^sphere — f 4 2) ’ 


■\y 


(15) 


The placement of the inclusion at x = 1/2 results in a local strain at the center of the inclusion 
equal to half the total, outer-radius model strain, which scales to half the model time in a 
simulation. 

We choose either a constant initial porosity field with (/)q = 0.05, or a random initial porosity 
with uniformly distributed values in the range 00±5 x 10“^. For the simulations with a randomly 
perturbed initial porosity field, we produced one initial field and reused this for all simulations. 
This initial field is created by first generating a random field on a uniform mesh that has a 
slightly larger grid size than the largest elements in the cylinder mesh; then this is interpolated 
onto the cylindrical mesh containing the spherical hole and variable grid size. This approach 
ensures that the random perturbations are sufficiently resolved by the mesh used in simulations 
and that the length scale of the perturbations does not vary with element size. 

Throughout this paper, we present simulation results on a two-dimensional slice through the 
inclusion at x = 1/2, as shown in Figure In this figure, the instantaneous compaction rate 
at time t = 0 for a simulation with bulk-to-shear-viscosity ratio i? = 20 is shown on the slice. 
The initial compaction rate is independent of the porosity exponent a for uniform porosity 
initial conditions. Pressure gradients caused by flow past the spherical inclusion induce two 
compacting lobes and two dilating lobes around the inclusion. This behavior was described 


in detail by McKenzie and Holness (2000); Alisic et al (2014); Rudge (2014) and is discussed 
further in Section [B. 11 


2.4 Discretization 

The problem described in Section |2.3| is solved by a finite element method on a mesh of tetra¬ 
hedral cells consisting of approximately 50 cells in the vertical dimension. The mesh is refined 
around the inclusion. The smallest cell size is ^ 3 x 10“^ near the inclusion, and the largest cell 
size is ^ 7 X 10“^ away from it. 

There are two main time stepping approaches to solving the two-phase flow equations 0 
namely, as a fully coupled system l\Katz et 'al\ |2007[) or by decoupling the porosity evolution 


equation 0 from the compaction equations ([2|)-([^ 
second approach. At each time step, equations 


{Katz and Takei, 2013). We follow the 


are solved to find the solid velocity, fluid 
pressure and compaction pressure, given the porosity and viscosities from the previous iteration. 
The porosity is then updated by solving 0 . To ensure a good approximation of the coupling, 
we iterate this process. Furthermore, if a non-Newtonian rheology is used, within each iteration 
a new strain rate is computed from the solid velocity and the viscosities are updated accordingly. 

The compaction system 0-0 is discretized with a continuous Galerkin finite element 
method using Taylor-Hood type elements (piecewise quadratic polynomial approximation for 
the solid velocity and piecewise linear polynomial approximation for the fluid and compaction 
pressures, see Rhebergen et al (2014)). The system of linear equations resulting from this dis¬ 


cretization is solved using Bi-CGSTAB in combination with the block-preconditioners developed 
in 


Rhebergen et al (2015). 


The porosity evolution equation 0 is discretized in space by a discontinuous Galerkin finite 
element method using a linear polynomial approximation. A Crank-Nicolson time stepping 
scheme is used to discretize in time (but using only the most recently computed velocity). 
To stabilize the simulation, a porosity-gradient-dependent artificial diffusion is added to the 
porosity evolution equation 0 of the form eV • (|V0pV0), with e = 0.1. To solve the resulting 
discrete system we use restarted GMRES preconditioned by algebraic multigrid. Simulations 
are terminated when the porosity becomes smaller than zero or larger than unity. 
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Instead of solving (§ @, it is possible to eliminate the compaction pressure by substitut¬ 
ing © into Q. The reduced system has fewer unknowns, but solving it is numerically less 
robust and less efficient than solving the expanded system We refer to \Rhehergen et aL\ 

(2015) for more details. 

Our simulation code is developed within the finite element software framework FEniCS/DOLFIN 


(Logg et a/., 2012 


solver library (Balay et al 


Logg and Wells 

2015a|bD . 


2010), in conjunction with the PETSc linear algebra and 


3 Results 

We group our results into two categories. In the first, the porosity is initially uniform. This 
means that the initial growth rate of the melt-banding instability is zero, and hence that changes 
in porosity are initially due solely to the presence of the inclusion. We consider the sensitivity of 
the compaction pattern to problem parameters, including the stress-dependence of the viscosity. 
The second category uses an initial condition with a random porosity perturbation. Melt-rich 
bands can potentially develop from the outset in this class of simulations. We explore in detail 
how and when such melt-rich bands develop, and how they interact with pressure shadows 
around the inclusion. 


3.1 Uniform initial porosity 


Newtonian viscosity We investigate the effect of the porosity exponent a and bulk-to-shear- 
viscosity ratio R on the porosity evolution in time-dependent simulations with a uniform initial 
porosity of 00 = 0.05 and Newtonian viscosity. When cr = 0 and n = 1, the shear viscosity is 
constant and uniform. In this case, the pressure shadows around the inclusion that are identified 
by perturbations in the porosity field rotate and advect with the matrix, with the top moving 
to the right and the bottom to the left, as shown in Figureand c (note that all cross-sections 
presented in this paper are oriented to have the same direction of shear). In contrast, in a 
simulation with a = 28 and n = 1, shown in Figure and d, the pressure shadows change 
shape in the opposite direction over time, following the orientation of expected bands in an 


inhomogeneous model (e.g. Spiegelman, 2003). 


To study the behavior of pressure shadows in more detail, we compute integrals of porosity 
on the two-dimensional slice through the inclusion. Integration is from the local radius of the 
edge of the inclusion r = a to one inclusion radius outward at r = 2a, for a series of azimuths 


between 0 and 27t around the circular cross-section of the inclusion (Qi et a/., 2013; Alisie et al 


2014) 


1 /• 2 “ 

“ / ' 

« Ja 


0dr. 


(16) 


Such radial integrals of porosity around the inclusion help expose the effect of a on the time 
evolution of pressure shadows. For the a = 0 simulation, the peaks become sharper over 
time, and the troughs become wider (see Figure ©)• The opposite happens for the cr = 28 
model (see Figure]^), with widening peaks. These differences are more pronounced for smaller 
bulk-to-shear-viscosity ratios i?, as seen in Figure These results are consistent with the 
two-dimensional results presented in Alisie et al. 


Non-Newtonian viscosity We introduce a non-Newtonian, power-law rheology in time- 
dependent simulations with uniform initial porosity. In these simulations, the power-law ex¬ 
ponent n is larger than one. The geometry of pressure shadows around the inclusion is affected 
by this strain-rate dependence, and ‘spokes’ form on either side of each pressure shadow quad¬ 
rant in a simulation with a = 0, as shown in Figure [^-b. This pattern is similar to the shape 


of pressure anomalies in non-Newtonian materials under simple shear found by Tenezer et al. 
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Figure 2: Results for simulations with a uniform initial porosity field and Newtonian shear 
viscosity (n = 1). The local strain at the center of the inclusion corresponds to one half of the 
reported model time, (a) Three-dimensional view of the porosity field for a simulation with 
porosity exponent cr = 0, bulk-to-shear-viscosity ratio i? = 5, at time t = 0.5 corresponding to a 
local strain of 0.25 at the center of the inclusion. The pressure shadows around the inclusion are 
shown as porosity contours of 0.045 in blue and 0.055 in red. (b) Three-dimensional view of the 
porosity field for a simulation with a = 28 and i? = 5, at time t — 0.5, with porosity contours 
at 0.045 and 0.055. (c) Slice through the porosity field at x = ^, for the simulation with a = 0, 
i? = 5, at t = 0.5. (d) Slice through the porosity field, for the simulation with a = 28, i? = 5, 
at t = 0.5. (e) Radial integrals over porosity, for the simulation with cr = 0, i? = 5, at various 
times t. (f) Radial integrals over porosity, for the simulation with a — 28, i? = 5, at various 
times t. 
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Figure 3: Radial integrals over porosity for simulations with uniform initial porosity field and n 
= 1, with various values of bulk-to-shear-viscosity ratio R. (a) Simulations with a = 0 at time 
t = 0.5, for various values of R. (b) Simulations with cr = 28 at time t = 0.25 or a local strain 
of 0.125 at the center of the inclusion, for various values of R. 


(2001). The second invariant of the strain-rate, shown in Figure]^ and which controls the vis¬ 
cosity variations, exhibits a complex pattern around the inclusion, without significant temporal 
variation throughout the simulation time. 

An increase in the power-law exponent n results in more pronounced spokes in the pressure 
shadows. However, these spokes mostly develop further than one inclusion radius away from the 
edge of the inclusion, and therefore the spoke shape is not reflected in the radial integrals (see 
Figure [^-e). Figure]^ further indicates that there is a decrease in amplitude of the peaks and 
troughs of the radial integrals for an increase in n. This implies that the strain-rate dependence 
of the viscosity does not enhance porosity growth rates. 

Increasing the porosity exponent a up to 28 (not shown here) does not result in a significant 
change of geometry of the pressure shadows in simulations with a total strain up to 0.2, indicating 
that the strain-rate dependence of the rheology is dominant over the porosity dependence at 
low strains. It is to be expected that at larger strains, when porosity anomalies have developed 
larger amplitudes, the porosity dependence becomes more significant. This could then lead to 
larger differences in geometry for a = 0 and 28. In contrast to porosity gradients, gradients in 
the strain-rate are large from the onset of simulations, as illustrated by Figure]^. The localized 
distribution of the strain-rate variations around the inclusion presents a significant resolution 
challenge for the numerical simulations. It has therefore proven difficult to model high strains 
for large values of n. 


3.2 Non-uniform initial porosity 

We now present simulations with a Newtonian rheology and initial porosity perturbations with 
a maximum amplitude of ±5 x 10“^ about a background porosity 00 of 0.05. The initial porosity 
field, shown in Figure [^, is the same for all simulations presented in this section. 

In a simulation with porosity exponent cr = 28 and bulk-to-shear-viscosity ratio i? = 1.7, 
melt-rich bands develop throughout the cylinder over time at an angle of ^ 45° with respect to 
the top and bottom of the domain (Figure [^). Larger band amplitudes are found towards the 
outside of the cylinder, as the local strain is proportional to the radius. Melt-rich bands develop 
both around the inclusion and away from it as shown in cross-sections at x = 1/2 through the 
inclusion and at x = —1/2 through the opposite side of the cylinder in Figure and c. In 
contrast, a simulation with the same a and i? = 20 does not display the formation of melt-rich 
bands, as shown in Figure and d. 

The integrals in Figure [^-f illustrate the difference in behavior between the R — 1.7 and 
i? = 20 simulations: the widening and flattening of the high-porosity peaks is much more visible 
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Figure 4: Results for simulations with a uniform initial porosity field and a non-Newtonian 
shear viscosity. The local strain at the center of the inclusion corresponds to one half of the 
reported model time, (a) Three-dimensional view of the porosity field for a simulation with 
porosity exponent a = 0, bulk-to-shear-viscosity ratio i? = 5, and power-law exponent n = 4 at 
time t = 0.5 corresponding to a local strain of 0.25 at the center of the inclusion. The pressure 
shadows around the inclusion are shown as porosity contours of 0.045 in blue and 0.055 in red. 
(b) Slice through the porosity field at x = ^, for the same simulation at t = 0.5. (c) Second 
invariant of the strain-rate field at t = 0. (d) Radial integrals over porosity, for the same 

simulation at various times, (e) Simulations with cr = 0 and i? = 5 at time t — 0.4, for various 
values of n. 
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Figure 5: Example of a simulation with a random initial porosity field, with a = 28 and R — 
1.7. (a) Initial porosity field, (b) Porosity field at t = 0.25. 
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Figure 6: Results for simulations with a random initial porosity field and Newtonian viscosity. 
The local strain at the center of the inclusion corresponds to one half of the reported model 
time, (a) Slice through the porosity field on the inclusion side of the cylinder at x = 1/2, for a 
simulation with a — 28, i? = 1.7, at t = 0.25. (b) Slice through porosity field on the inclusion 
side of the cylinder, for a simulation with a = 28, i? = 20, at t = 1.0. (c) Slice through 
the porosity field on the side of the cylinder opposite the inclusion at x = — with a — 28, 
R = 1.7, at t = 0.25. (d) Slice through the porosity field on the side of the cylinder opposite the 
inclusion, with a = 28, i? = 20, at t = 1.0. (e) Radial integrals over porosity, for the simulation 
with a = 28, R — 1.7, at various times t. The solid lines are fits with Fourier functions with 
the lowest nine coefficients included, (f) Radial integrals over porosity, for the simulation with 
cr = 28, i? = 20, at various times t. 
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Figure 7: Radial integrals over porosity for simulations with a random initial porosity field 
and Newtonian viscosity, with various values of bulk-to-shear-viscosity ratio R. The solid lines 
are fits with Fourier functions with the lowest nine coefficients included, (a) Simulations with 
cr = 15 at t = 0.5, for various values of R. (b) Simulations with a = 28 at t = 0.5, for various 
values of R. 


in the i? = 1.7 case than in the i? = 20 case. In the latter case, the porosity shadows even display 
an advected pattern at large strains (represented by sharp peaks much like the simulation with a 
uniform initial porosity field in Figure]^), indicating that the growth of porosity is less dominant 
than its advection for such large R. 

Melt-rich bands only develop in simulations with sufficiently large a and small i?, as illus¬ 
trated by the more pronounced widening and flattening of high-porosity peaks in the integrals 
in Figure This is in line with the expected growth rates of melt-rich bands derived using 
linear stability analysis and presented in Appendix |B.2[ The linear stability analysis predicts 
melt bands to grow initially exponentially (oc exp{st)) at a dimensionless rate 


s = 


0^(1 - 0 q ) 


(17) 


which indicates that melt-rich bands are expected to grow faster for larger a and smaller R. 


3.3 Model regimes 

Figure [^summarizes the results of our parameter study of porosity exponent a and bulk-to-shear- 
viscosity ratio R for simulations with a random initial porosity field and Newtonian rheology. 


The overall pattern is similar to that found in the two-dimensional study of Alisic et al. (2014): 
melt-rich bands only develop for i? < 5 and large a. The region of the parameter space in 
which bands develop is slightly larger for the three-dimensional geometry compared to a two¬ 


dimensional case (see Alisic et al (2014, Figure 10)), and the maximum strain achieved in the 
simulations is significantly larger. This might be explained by the fact that the amplitudes of 
pressure shadows decay faster away from the inclusion in three dimensions compared to two 
dimensions (as rather than Rudge (2014)), resulting in less dominant pressure shadows 


compared to other features developing in the porosity field. 


4 Discussion 


The simulations in this study display several of the main features observed in the laboratory 
experiments by Qi et al. (2013), such as the pressure shadows around the spherical inclusion. For 
small bulk-to-shear-viscosity ratios, melt-rich bands develop throughout the medium, including 
in the vicinity of the inclusion. However, these bands do not have the dominance observed in 
the laboratory experiments, where they grow as very straight and pervasive features directly 
adjacent to the inclusion, overprinting the pressure shadows around the inclusion. In contrast. 
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Figure 8: Summary of parameter study of bulk-to-shear-viscosity ratio R versus porosity ex¬ 
ponent a. Red circles indicate that no significant development of melt-rich bands takes place 
during a simulation with the specified combination of R and a; green squares indicate the pres¬ 
ence of melt-rich bands. Blue triangles indicate development of melt-rich bands only away from 
the inclusion. The black contours denote the maximum strain achieved at the outer edge of the 
cylinder before the porosity goes out of bounds and the simulations end (0 < 0 or 0 > 1). 


the bands found in the numerical simulations grow to shorter lengths and do not overprint the 
porosity structure around the inclusion. Furthermore, simulations with and without random 
initial porosity perturbations using a given combination of a and R and a Newtonian rheology 
run to the same maximum model strain before going out of bounds. This behavior, along with 
the fact that melt bands form for a larger (cr, R) parameter space away from the inclusion 
than near it, are indications that in our models the pressure shadows around the inclusion are 
dominant over any bands that form in simulations with random heterogeneities. Several recent 


studies investigate alternative constitutive relations (Takei and Katz, 2013; Katz and Takei 


2013; Rudge and Bercovici, 2015) that could potentially affect the balance of pressure shadows 


and melt-rich band formation near the inclusion. 

The dominance of pressure shadows also points to a key deficiency in current models of two- 
phase-flow: the models contain no physics that sufficiently limit porosity growth, which results 
in the porosity field in our models growing until reaching unity, at which time the simulations 
are terminated. Realistically, the governing equations are only valid for porosities much smaller 
than unity. A second consequence of the porosity weakening rheology in our model is the lack 
of a minimum length scale (width) to which melt-rich bands will evolve. This means that the 
thickness of bands in simulations is ultimately dictated by the grid spacing, and therefore the 
solutions are resolution dependent. 

It should be noted that our numerical simulations have a different velocity boundary condi¬ 
tion on the sides of the cylinder than the experiments by Qi et al (2013). In those experiments, 
velocity is only prescribed on the top and bottom of the cylinder, and the side boundary can slip 
freely. In contrast, in our simulations the velocity is fully prescribed on all outside boundaries, 
leading to a potentially more constrained model. 

In simulations with a rheology that is also strain-rate-dependent, limitations on numerical 
resolution near the inclusion prevented evolution to high strains. Therefore the regime where 
amplitudes of porosity variations were large enough to allow the porosity-weakening to become 
dominant over strain-rate effects was typically not reached, and bands would not develop within 
the simulation time. 

In future numerical studies, it would be helpful to utilize a significantly higher resolution 
near the inclusion in such simulations, so that the strain-rate gradients can be better resolved. 
In addition, much could be gained from obtaining higher-resolution images of pressure shadows 
in experiments: the details of the shape of the shadows could help with identification of the 
prevailing deformation mechanism (diffusion creep or dislocation creep with large power-law 
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exponent n). 


5 Conclusions 

We have modeled the behavior of partially molten material with an inclusion under torsion using 
three-dimensional numerical solutions to the equations of two-phase flow. Recent laboratory 
experiments with a similar setup display a competition between pressure shadows forming around 
the inclusion and melt-rich bands that develop throughout the partially molten medium. In our 
numerical simulations, the pressure shadows around the inclusion are reproduced for all tested 
combinations of bulk-to-shear-viscosity ratio and porosity exponent of the shear viscosity. In 
contrast, melt-rich bands only develop for small bulk-to-shear-viscosity ratios of five or less. We 
conclude that it is more difficult to form melt-rich bands near the inclusion, which provides a 
strong perturbation to the pressure field in the form of pressure shadows. Comparing this study 
with our earlier work in two dimensions, we show that the pressure shadows are less dominant 
in three dimensions, resulting in more pervasive development of melt-rich bands. For strain-rate 
dependent viscosity, the shape of the pressure shadows is significantly different compared to 
Newtonian cases. This variation in shape could be utilized to pinpoint the dominant deformation 
mechanism around the inclusion in future experiments. 


A Pressure splitting and nondimensionalization of the govern¬ 
ing equations 


The dimensional equations for two-phase flow are: 

^_V.(l-0)u, = O, 

V • u = 0, 

(p{uf - Us) = -^Vpf, 

M/ 

V O’ = 0, 


(18) 

(19) 

( 20 ) 
( 21 ) 

where 0 denotes porosity, t time, and and uj the solid and fluid velocities, respectively. Bulk 
properties are denoted with an overbar, where a bulk quantity a = (/)aj-h(l —(/))a 5 . Furthermore, 
Kfj) is the permeability, /ij the fluid viscosity, pf the fluid pressure, and a is the bulk stress tensor. 

We define the bulk stress tensor in terms of the fluid pressure pj, compaction pressure pc, 
and the deviatoric stress tensor r: 


a = -pfl-pcl + T, 

Pc = -CV • u^. 


r = 2pe = p Vu^ -h (Vu^)^ - - (V • u^) I 


( 22 ) 

(23) 

(24) 


where I is the identity tensor, ( the bulk viscosity, p the shear viscosity, and e the deviatoric 
strain rate tensor. 

We can now write a new system of equations using u^, pj, pc, and 0 as unknowns: 


^ - V • (1 - (/>) u, = 0, 

(25) 

+ 

< 

< 

II 

(26) 

\ / 

-V • Us - C“Vc = 0, 

(27) 

-V • T + Vp/ + Vpc = 0- 

(28) 
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Constitutive properties are defined in this study as follows: 


^(f) — -^ref 



V = ^ref 



■q;(0-(/)o) 

5 


C — Cref 



a{<p—<po) 




(29) 

(30) 

(31) 


where n = 2 , m = 1, and a is the porosity exponent and denotes the permeability at the 
reference porosity ^o- The second invariant of the deviatoric strain-rate tensor i is given by: 

/I 

’ ( 22 ) 

and q is related to the power law exponent n by 

g = l-t. (33) 

The reference value of the second invariant is chosen as 

. _ ip 

“ 2H 

which is the value the second invariant takes on the curved boundary of the cylinder under the 
imposed torsion field. H is the radius of the cylinder and 7 is the imposed shear strain rate on 
the curved boundary, r/ref smd thus represent the shear and bulk viscosities at the curved 
boundary when the porosity is uniform and equal to 00 . The bulk-to-shear-viscosity ratio R is 
given by Cref/href- 

To complete the problem, boundary conditions are applied as follows: 

- -'Vpf • n = 0 on (35) 

fif 

U 5 = w on (36) 


(34) 


where w is a prescribed solid velocity, and the boundaries are taken to be impermeable. 

We must now define a convention for nondimensionalizing the governing equations, using 
primes for dimensionless quantities: 


X = Hx.', Us = HjUs, t = j 

Pf = rireap'f, Pc = hrefTTc, (37) 

^<f> h hrefh i C CrefC • 

The dimensionless system of equations then becomes, after dropping the primes, equations 
0-0 in the main text. 


B Analysis and code benchmarks 

B.l Instantaneous compaction around an inclusion 

The instantaneous rate of compaction around a spherical inclusion in an unbounded medium 


with uniform porosity can be derived analytically [McKenzie and Holness, 2000; Rudge, 2014), 


which therefore lends itself for benchmarking our numerical method. 
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In Rudge (2014) analytical solutions were presented for compacting flow past a sphere in 


far-held simple shear. Here we generalize this solution to the case of a far-held torsional how. 


We follow the approach of Rudge (2014) in this appendix, which has a different set of coordinates 


to the main body of this paper, with the center of the sphere being the origin of the coordinate 
system. In Cartesian coordinates the far-held torsional how takes the form 


ur 


= r {-{y - yo){z - zo), (x - xo){z - 2 :o), 0 ) , 


(38) 


where the origin of the torsional how is at (xq, 2/o, ^o) sind F is the twist rate. The above can be 
decomposed into irreducible Cartesian tensors as 




— F2 X X “h E 


1 


X — 


-X 


{6 


X) 


where 


V = r {-yoZQ, xqzq, 0 ), 


(39) 

(40) 

(41) 

(42) 

(43) 


and X = {x, y, z) is the position vector. According to the Faxen laws a sphere placed at the 
origin in such a how will translate with velocity V and rotate with angular velocity i7, provided 
that there is no net force or torque on the sphere. The compacting how past the sphere can be 
calculated as a linear superposition of the how due to pure strain E • x ^Rudge , 2014, Sec tion 5) 
and that due to the vortlet how —lx x {6 • x) (a quadratic how, not considered in Rudge 



( 

XQ yo 

2 ’ 2 ’ ' 

-^oj 

5 



( ° 


0 

yQl‘^ 

\ 

r 

0 


0 

-xo/2 



yo/2 

-X0I2 

0 

/ 



/ 

-1 

0 

0 ^ 



0: 

= r 


0 

-1 

0 
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\ 

0 
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(2014)). The vortlet how is characterized by the second-rank pseudo-tensor 6 which is equal to 
the vorticity gradient. The perturbation how satishes boundary conditions 


Us|r=a, — 2 ^ ^ ' 

Us ^ 0 , 


X 


^Pf • nl^.—Q, — 0, 


(44) 

(45) 


pj ^ 0, as r ^ oo, 

where r = \x\ is distance from the center of the sphere, and a is the radius of the sphere. The 
solution to the governing equations with these boundary conditions does not involve compaction 
and is simply a Stokes how, given by 

U3 = ^x X (e-x), (46) 

Pf = 0. (47) 

Since the quadratic how does not involve compaction, the instantaneous compaction rate 
for a sphere in a torsional held is the same as that for pure and simple shear, given by Rudge\ 
(2014, eqns. (5.51) and (5.63)). When the compaction length is large compared to the domain 
size, the behavior can be well described by the large-compaction-length asymptotic limit of the 
equations, where the instantaneous compaction rate and huid pressure are given by 

V-U, = (48) 


Pf 


_ 


2z/ “h 3 
5z/ 


i^ref 6 (2z/ -\- 3) 


3a 


X • E • X 


(49) 
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Figure 9: Instantaneous pressure fields for a simulation with inclusion size 0.1, porosity exponent 
a = 28, bulk-to-shear-viscosity ratio i? = 5/3, and stress exponent n = 1. (a) Compaction 
pressure. (b) Fluid pressure. 


where u = r^ref/ (Cref + 477r ef/3) (i^ + 4/3) Th e above solution is identical to that given in 

equations (30) and (32) oi \McKenzie and Holness\l\2000l . 

In Figure]^ we show the numerically-calculated instantaneous compaction pressure and fluid 
pressure on the two-dimensional slice indicated in Figure [^, resulting from the imposed torsion. 
We use a large compaction length with D = 100, and the same cylindrical mesh as outlined in 
Section [23| The numerically-calculated compaction pressure, shown in Figure]^, very closely 
matches the analytical expression in (48). The fluid pressure, shown in Figure]^, matches less 
well, with the expected quadrupole pattern of (49) disturbed at the top and bottom boundaries 
of the finite computational domain. That the fluid pressure is affected more by the boundaries 
than the compaction pressure is to be expected from the analytical expressions in (48) and (49); 
compaction pressure decays rapidly away from the inclusion, as whereas fluid pressure 
decays much more slowly, as and thus the fluid pressure feels effects at larger distances. 

For further validation of the numerical simulations, we compute L 2 error norms for the 
fluid pressure pf^ compaction pressure pc, and solid velocity u^, with respect to the analytical 
solutions. We define the following error for a field y: 




\\x^-x' 

llx^ll 


(50) 


where the numerical field is denoted by and the analytical solution by We compute 
this for a series of inclusion radii between 0.05 and 0.2, as shown in Figure [T^ The L 2 error 
decreases with decreasing inclusion radius. This indicates that the error with respect to the 
analytical solution results from the presence of boundaries in the numerical domain, which do 
not exist in the analytical solution. This effect becomes less dominant when the inclusion is 
further away from the outside cylinder boundaries, i.e., for smaller inclusions. In addition, the 
L 2 error is larger for the fluid pressure than for the compaction pressure. This is due to the 
same boundary effects as observed in Figure 


B.2 Linear stability analysis of melt bands under torsion 

Linear stability analysis provides important insight into the expected growth rate of melt bands 


(e.g. 

Stevenson, 1989 

; Spiegelman 

, 2003; 

Katz et al, 2006; 

Butler, 2009 

; Takei and Katz 

, 2013; 

Rudge and Bercovici 

, |2015). The solutions arising from linear stability analysis can also be 


used as a check on numerical solutions of the full set of governing equations (see Alisic et al 


(2014, Appendix C2)). In this appendix we present a linear stability analysis of melt bands 
under torsion in an infinite cylinder for our chosen rheology, and use the solutions to benchmark 
our numerical code. The analysis closely follows an earlier linear stability analysis of melt 


bands under torsion by Takei and Katz (2013). In Takei and Katz (2013) it was assumed that 
perturbation wavenumbers are large, which allows the neglect of various radial derivatives in 
the analysis. We do not make this assumption here, and instead solve numerically for the radial 
variation. 
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Figure 10: L 2 error norms computed for simulations with different inclusion radii, for fluid 
pressure pf^ compaction pressure pc, and solid velocity perturbation Aug with respect to a 
torsional velocity field without an inclusion. 


The base state solution has a uniform porosity 00 and solid velocity uq given in cylindrical 
coordinates (p, 0, z) as 

Uq = (51) 

where F is the twist rate. The twist rate F is related to the shear strain rate 7 on the cylinder 
edge by 7 = The base state solution has zero pressure everywhere (po = 0), is not com¬ 

pacting such that Co = V • uo = 0, and has a strain-rate-tensor with only the (0, z) component 
non-zero, 

eo = r^(^^z +. (52) 

We seek small perturbations about this base state of the form 


^ = 00 + e(/>i + • 
Us = uo + eui + 
Pf = 0 + epi + ■ ■ 


(53) 

(54) 

(55) 


Substituting (53)-(55) into the governing equations (18)-(21) leads to equations at first order 
in e given by 


= (1 - 0o)Ci, (56) 

Cl - — VV = 0, (57) 

M/ 

V • 0-1 = 0, (58) 

0-1 =-pil + (oCiI + 2rjoei + 2pieo, (59) 


where ^ | + uq • V. 

When the rheology is non-Newtonian, the base state viscosities vary with radius according 
to 


Vo = Vr:ef {P/H) , 

C0 = Cref 


Expanding the rheological law (30) to first order in e yields 


m = Vo 


^-o;0i - q 


ei : eo\ 
) 




(60) 

(61) 


(62) 
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The expected growth rate of bands can be determined by replacing (56) by 


s4>i = (1 - </>o)Ci, 


( 63 ) 

and ([6^-([63|) for the instantaneous 


and solving the eigenvalue problem described by 
growth rate s. The finite element method can be used to numerically solve for the eigenfunctions 
and eigenvalues of the linear stability equations. Equations ([^-([^ and (62)-(63) can be cast 
into a weak form for trial functions {u^p) and test functions {v^q) as 

2r?oe“ : e" + CoCV" - pC^ - qC^ - yVp-Vq - dV = A ^ 4r?oeoC“e;, dF, 


e«^i(^Vu + (Vu)^)-4(V-u)I, 

= V • u, 


(64) 

(65) 

( 66 ) 


where the subscripts 1 referring to the first order state have been neglected for clarity. The 
eigenvalue A is related to the growth rate s by 


A = 


a(l - ^o) 


(67) 


The impermeable and no-slip boundary conditions (35) and (36) on the cylinder edge lead to 


the vanishing of surface integral terms in the weak form (64). 


The three-dimensional weak form (64) can be reduced to a weak form for the radial coor¬ 
dinate alone using symmetry considerations. Invariance of the cylinder under rotation about 
its axis, and invariance under translation in the 2; direction, suggests looking for eigenfunctions 
proportional to where n is the angular wavenumber, and h is the vertical wavenumber. 

Substituting solutions of this form 



(68) 


(69) 


(70) 


(71) 


and corresponding complex conjugates for test functions into (64) leads to purely radial integrals 
where dV 27rpdp. The integrands are real, and the resulting radial eigenfunction problem 


was solved using FEniCS/DOLFIN (Logg and Wells, 2010| and the eigenvalue solver SLEPc 


{ Hernandez et a/.]|2005| ). An example eigenfunction calculated using this approach is shown in 
Figure 11 


B.2.1 A special case: Newtonian melt bands under torsion in an infinite domain 

There is a special case of the linear stability analysis for which a complete analytical solution 
can be obtained. This is the case when the rheology is Newtonian (n = 1), and the domain 
of interest is infinite. In this case the solenoidal component of the flow is decoupled from the 


irrotational component, simplifying the analysis (Spiegelman 
equations are 


1993, 2003). The linear stability 


dt ^ 


(1 - 4>q)C\, 


-V^Ci + S-^Ci = -2uar^^, 

ozoyj 


(72) 

(73) 
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Figure 11: Porosity field for a particular eigenfunction of the linear stability equations. This 
eigenfunction has angular wavenumber n = 5, and vertical wavenumber h chosen such that the 
angle of the bands to the shear plane is 45° on the edge of the cylinder. The bulk-to-shear- 
viscosity ratio i? = 5/3 and the power law exponent n = 1. The compaction length is large 
(100 times the cylinder radius). The eigenfunction shown is the fastest growing mode with this 
choice of n and h. 


where (72) follows from (56), and (73) is a result of combining the divergence of (58) with (57), 
(59), and (62), 6 is the compaction length and u = T^ref/ (Cref + 477ref/3). Solutions to ( |72[ ) 
and (73) can be found in the form of cylindrical harmonics (eigenfunctions of the Laplacian 
operator in cylindricals), 


(74) 

Cl = (75) 

where Jn{z) is a Bessel function of the first kind and h{t) is the vertical wavenumber, which 
varies with time as 

h{t) = /lo — Frit^ (76) 

due to the advection. Note that 


V^Ci = -k'^{t)Ci, 
k^{t) = \^ + h^{t). 


Substituting (74) and (75) into (72) and ( |73[ ) yields 

m = (1 - (i>o)c{t), 

S~‘^ + C(t) = 2i/arnh(t)^(t), 


which can be combined to give 




(77) 

(78) 

(79) 

(80) 


(81) 


and integrated to give 

(( 5“2 + J 



(82) 


These expressions closely mirror the expressions for simple shear given by 
compare (81) and (82) here with equations (27) and (33), respectively, from 


Spiegelman 

(2003) 

Spiegelman 

(2003) 
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Figure 12: L 2 error norms of growth rate at the center-plane of the cylinder, computed for 
numerical simulations with different resolutions. The eigenfunction tested is shown in Figure 

ttH 


The expressions are identical, except with the appropriate switch of angular wavenumbers for 
planar wavenumbers, and the difference in nondimensionalization. 

In a cylinder of finite radius, common choices of boundary conditions on the cylinder edge 
(including no-slip) lead to a coupling of the solenoidal component of the ffow and the irrotational 
component, which complicates the analysis given above. Nevertheless, this special solution can 
be used as a check on the numerical approach to calculating the eigenfunctions described in the 
preceding section. The special solution (74) and (81) was recovered numerically for a choice 
of boundary conditions that decouples the solenoidal ffow from the irrotational. This choice 
of boundary condition is impermeable, and almost, but not quite, free-slip: on the cylinder 
edge Up — 0, Gpz — 0, ap^ip = —2z/?i^/p, and dpf/dp = 0. These boundary conditions imply 
that dC/dp = 0 on the cylinder edge, and restrict the allowable values of A to the roots of the 
derivative of the Bessel function. As is the case for simple shear, the fastest growing modes 
occur for infinite wavenumbers. The maximum growth rate occurs as n ^ oo, h n/H (i.e., 
bands at 45° to the shear plane on the cylinder edge) where s vol{\ — (/)o)Ti7. 


B.2.2 Linear stability benchmark 

We test the application code by numerically computing the instantaneous growth rate of porosity 
for an initial porosity field given by the eigenfunction shown in Figure The numerically 
computed growth rate can be compared with the expected growth rate of the eigenfunction 
determined from the eigenvalue. 

Figure [T^ shows an example of such a comparison, where an error norm for growth rate is 
plotted against resolution for various mesh resolutions from approximately 10 to 50 elements in 
the vertical direction, to study the effect of grid size on the accuracy of the numerical method. 
It is important to note that the eigenfunctions are determined for a cylinder of infinite extent, 
whereas the simulation domain is a cylinder of finite extent. To mitigate the resulting boundary 
effects, the two are compared only on a slice through the center-plane of the cylinder, at z = 1/2. 
The L2 error norm is calculated for the local instantaneous growth rate of porosity on the slice, 
using equation (50). Generally, the error in growth rate on the center-plane decreases with 
increasing resolution (i.e., with decreasing grid size), until a limit is reached at a grid size 
around 0.03. For finer grids the error does not decrease any further, which we attribute to the 
effect of the top and bottom boundaries on the growth of porosity. Computed growth rates are 
typically within a few per cent of the expected growth rates, which gives us confidence that the 
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application code is solving the compaction equations effectively. 
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